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Abstract 

Crack propagation is studied numerically using a continuum phase-field approach to mode III 
brittle fracture. The results shed light on the physics that controls the speed of accelerating cracks 
and the characteristic branching instability at a fraction of the wave speed. 
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The quest for a fundamental understanding of dynamic brittle fracture has been an 
ongoing challenge. The traditional continuum approach to computing the speed of brittle 
cracks consists of solving the equations of linear^lasticity with boundary conditions on 
the moving fracture surfaces up to the crack tip [1]. The solutions have stress fields that 
diverge near the tip representing a finite energy flow rate to the tip. The crack speed v is 
then assumed to be uniquely determined by this energy flow rate. In this theory, all the 
nonlinear physics of failure inside a microscopic region around the tip — the so-called process 
zone — is buried in a phenomenological function that relates the fracture energy F of the 
material (energy needed to advance the tip per unit length of crack front and per unit of 
crack extension) and v. 

Detailed experiments have been carried out to test this theory by measuring the speed of 
accelerating cracks with different initial lengths and different loads Q]. The theory should be 
valid if a unique T{v) curve applies to the different cracks for the same material. The main 
difficulty in any experimental test of this theory is that the energy flux to the tip cannot be 
measured directly. Therefore, this flux needs to be inferred from a time-dependent solution of 
linear elasticity for accelerating cracks. Using an approximate solution, Sharon and Fineberg 
have concluded that the data for different cracks collapse on a unique F(f) curve, whereas, 
using an exact solution of Eshelby [3j , Kessler and Levine have concluded that this collapse 
does not occur [J]. 

Also of long-standing interest in brittle fracture is the existence of a dynamic instability 
that limits the speed of fast moving cracks. This generic instability has been seen in both 
expeH^en. and .olecla. ,y.^,. ..ula.io. Q. The tac. ..at d.ffe.n. a„„. 

phous materials with entirely different microscopic details (such as glass and PMMA) exhibit 
strikingly similar branching instabilities S] strongly suggests that a continuum theory may 
be appropriate for understanding this phenomenon. Devising such a theory, however, has 
proven to be difficult. Cohesive zone theories modify the boundary conditions on the stress 
field near the tip to take into account the short scale force between crack surfaces. Langer 



L„b.ovs.y B, .oweve. Kave sKown that tKese ,node>s a. unsuitable f„. stab„,t. cal- 

culations since the results depend singularly on the details of the cohesive zone. In addition, 
intrinsically discrete branching instabilities in lattice models 0, seem unlikely to bear 
relevance to experiments in amorphous materials. 

In this letter, we study dynamic brittle fracture using a recently developed continuum 
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phase-field approach for mode III cracks 



This approach has the chief advantage that it 
incorporates both the short-scale physics of failure and macroscopic linear elasticity within 
a self-consistent set of coupled nonlinear partial differential equations that can be solved 
numerically. Moreover, solutions free of discretization artifacts can be obtained. Earlier 
simulations of this model were limited to small system sizes owing to the extreme stiffness of 
the equations, and no branching instability was found even for cracks approaching the wave 



speed ^|. Here, we briefly discuss how to overcome this stiffness and carry out simulations 
in systems large enough to make contact with results from the fracture community. In 
addition, these simulations enable us to study the onset of the characteristic branching 
instability at a fraction of the wave speed that turns out to be a robust feature of this 
phase-field model. 

The basic va.ab.es of *e ,„odel Q a.e the scala. disp.aee.e„. „ pe.pendieu.a. 

to the x-y plane of mass points from their original positions, and the phase-field, (f){x,y), 
which describes the state of the material. The unbroken solid, which behaves purely elas- 
tically, corresponds to = 1, whereas the fully broken material that cannot support stress 
corresponds to = 0. The total energy (kinetic plus elastic) of the system per unit length 
of the crack front is 
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E = j dxdy 



f (I) +f |V0p + /^/(0) + f^?(0)(|ef -e^) 



(1) 



where p is the density, t = Vm is the strain, /(0) is a double- well potential with minima at 
= 1 and = 0, /i is the elastic shear modulus, and ec is the critical strain magnitude such 
that the unbroken (broken) state is energetically favored for |e| < Cc (|e| > ec). The function 
g{(f)) is a monotonously increasing function of with limits g{0) = and ^'(0) = 1, which 
controls the softening of the elastic energy as more bonds are broken. 

Taking the first variations of the energy with respect to the strain and to 0, we obtain 
the stress, a = 6E/6e, and the equations of motion 

90 5E 
9^ = 50' 

P|^-V... (2b) 

Energy is dissipated in the process zone around the crack tip where varies rapidly in 
space and time. Eq. ^la^ implies that the size of the process zone is ~ ^ = k/ pel and 



the characteristic time of energy dissipation in this zone is r = l/(x/ie^). Furthermore, by 
rescahng lengths by ^, time by ^/c, where c = a//Vp is the shear wave speed, and m by ,^ec 
in Eqs. we find that crack propagation is controlled by two dimensionless parameters 



6 = h/{fiel) and /? = cr/^. The first determines the surface energy normalized by fie^C, jll] 

^ /' Vl-^(</.) + 25/(0). (3) 
Jo 
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For the choices /(0) = 160^(1 — 0)^ and g{(f)) = 40^ — 30^ here, 7 approximately doubles 
when 6 increases from to 2. 

Parameter (3 controls the importance of inertia relative to dissipation in the process zone. 
When /3 -C 1 (inertia-dominated regime), failure is rapid and crack propagation is governed 
by elastic energy fiow. When j3 ^ 1 (dissipation-dominated regime), failure is sluggish. 
The elastic displacements are quasistatic and propagation is governed by the kinetics of the 
failure process. 

We study fracture numerically in a strip of width 2W with a fixed displacement 
u{x, ±W) = ±A at the strip edges. The stored elastic energy per unit area ahead of 
the crack tip is G = fiA'^/W and the Griffith's threshold load for a semi- infinite crack is 
Gc = 27. 

The equations of the model are inherently stiff due to strain localization behind the crack 
tip The strain is localized within a narrow region whose width vanishes in the large 
strip limit W/^ 00, while the width of the variation remains of order ^ in this limit. 
Therefore, while it is possible to resolve the spatial variation of everywhere, resolving 
the strain concentration behind the crack tip is computationally impractical for large strip 
widths. The key to overcoming this difficulty is to recognize that the model remains well- 
defined when the strain variation is resolved only in the process zone, but is allowed to be 
discontinuous on the lattice scale far behind the tip. Therefore, it is possible to simply 
discretize the energy E on a square lattice and to formulate the dynamics of the model from 
partial derivatives of E with respect to the variables and u at lattice points. One subtlety 
of the present simulations is that high-frequency waves are radiated behind the crack tip 
above a small threshold tip speed. These waves increase slightly the fracture toughness and 
originate from the region behind the tip where the displacement becomes discontinuous on 
the lattice scale. In addition, multiple refiections of these waves can cause tip oscillations 
at long times. This time, however, is longer than the time for propagation to reach a 
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FIG. 1: Crack velocity vs. tip position for cracks accelerating from rest with different initial lengths. 
P = 1, 6 = 0, lattice spacing dx = 0.3^, and W = 30^. The velocity for the largest load becomes 
negative since the symmetric branches retreat to give way to the asymmetric tip-splitting mode 
with a sinusoidal (snake-like) fracture path. 

steady-state or branching to occur such that this feature is not problematic. Furthermore, 
this radiation is easily suppressed by adding a small Kelvin viscosity to Eq. ()2bj) with no 
qualitative change of the results presented here. 

We have verified that the important observables converge quickly in the limit of the 
vanishing lattice spacing, and that we operate in the regime where these observables differ 
by at most 15% from their continuum limits for the lattice spacing dx in the range 0.3^ to 
0.4e 

To study unsteady crack propagation, we use stationary solutions of Eqs. Q as initial 
conditions that correspond to cracks of different initial length at rest. For a given load G 
and corresponding initial crack length i, these stationary solutions are found by relaxing u 
with the Gauss-Seidel iteration scheme, which amounts to solving Eq. ()2b|) without inertia, 
while is relaxed using Eq. (j^aj) . This procedure yields stationary cracks whose initial 
length i decreases with increasing G. We then simulate the full equations of motion with 
inertia to study accelerating cracks with zero initial velocity. These simulations parallels 
previous experiments [3] and lattice simulations ^\. To study steady-state features of crack 
propagation independent of initial conditions, we run long simulations in strips that are 
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FIG. 2: Fracture energy vs. velocity for the accelerating cracks and the same parameters as in Fig. 
n The solid line corresponds to steady-state crack propagation with F = G by energy conservation. 
The five = 1/2 contours shown on the inset correspond to the five large solid triangles. They are 
equidistant in time separated by 10^/c. Blunting and subsequent splitting is responsible for the 
deceleration. 

effectively infinite along the propagation direction. To keep the computations tractable, 
we periodically translate the fields u and (p by one lattice spacing such that the crack tip 
remains in the middle of a strip of length much larger than W. We have checked that the 
results of these "tread mill" simulations are independent of boundary effects. 

To compute the fracture energy during unsteady motion, we equate vT with the expression 
for the energy flow rate to the tip of mode III cracks jsl 

vr= dC[fiudnU + vn^dt{puy2 + fi\e]'^/2)], (4) 
Jc 

where C is a closed circuit around the moving tip and n is the outward normal to C. Since 
energy is conserved in the model where = 1, we can equivalently obtain the energy flow 
rate to the tip by calculating the time rate of change of the total energy energy (elastic plus 
kinetic) in the region of the system where is larger than a threshold value 0c arbitrarily 
close to unity. This quantity is precisely the energy flow rate into the process zone defined 
as the region where < 0c. 

Plots of tip speed versus tip position for cracks accelerating from rest are shown in Fig. ^ 
for the inertia-dominated regime (3 = 1. Corresponding plots of F versus v are shown in 
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FIG. 3: Contours of = 1/2 separated in time by At = 10.^/c. All plots are drawn on the same 
scale for the parameters P = 2, 6 = 0, and W = 30.^. (a) Transient branching for G = 1.56Gc. (b) 
Weak periodic branching (the "snake") for G = 1.86Gc. (c) Chaotic branching for G = 2.90Gc. 

Fig. El The initial crack acceleration increases with load and the crack tip splits into two 
symmetric branches above a critical onset load Gonset- Cracks for smaller loads which do not 
split, or split only transiently for G ~ Gonset, reach a steady-state propagation velocity that 
coincides with the steady-state velocity calculated on the tread mill (solid line in Fig. 
The fact that F = G in steady-state is a self-consistency check of our method of calculating 
F for unsteady cracks. 

From the long simulations on the tread-mill, three basic regimes of crack propagation can 
be distinguished: a stable regime for small loads, where v is constant in time and the crack 
is rectilinear, an asymmetric tip-splitting ("snake") mode for intermediate loads, where v 
fluctuates in time around some average value while the crack follows a sinusoidal trajectory, 
and a chaotic tip-splitting regime for large loads. Examples of these regimes are shown in 
Fig. ini The simulations also reveal a strong coupling between the crack tip dynamics and 
elastic waves that we have not explored in detail. A clear signature of this coupling is the 
fact that the temporal period of branching in the snake mode is very close to the period 
AW/c of the first harmonic standing wave of the strip. The transient coupling of the tip to 
higher frequency waves is also seen for higher loads in the chaotic branching regime. 

The steady state velocity v saturates at some value Vc- Rectilinear cracks cannot propa- 
gate faster than Vc- We conclude that off-axis branching in the present model is due to the 
absence of steady-state crack solutions above a critical speed, rather than a linear instability 
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FIG. 4: Limiting speed of crack propagation Vc vs. f3 for different values of 6 (or 7 given by Eq. 
131). Inset: in the inertia-dominated (/? — > 0) limit vs. 5. 

of solutions that exist up to the wave speed, as in the lattice models [l^. The dependence 
of Vc on the parameters of the model is shown in Fig. EJ The maximum crack speed is 
well defined in the inertia dominated limit (3^0. It grows monotonically with the scaled 
surface energy 7 (Eq. Q) and has a minimum at 6 = oi Vc ^ 0.41c. At this speed the 
linear elastic field in the unbroken material around the crack tip is quasi-isotropic We 
therefore conjecture, along the lines of Gao jisl, that tip blunting which leads to velocity 
saturation and ultimately to tip-splitting is due to the relativistic contraction of stress fields 
in the nonlinear process zone where the sound speed is small due to the softening of the 
material. 

The smallest value of Vc ~ 0.41c in the (3 —>■ limit is consistent with the calculation of 
Adda Bedia [14J which shows that tip splitting is energetically possible for mode III fracture 
for speeds above 0.39c with a maximum angle between symmetric branches of 80 degrees, 
which is about 10-15 percent larger than the angle in our simulations. It should be noted, 
however, that this calculation 14] considers a steady-state crack that stops abruptly and 
splits into two daughter cracks, whereas cracks decelerate before splitting in our simulations. 
In addition, the onset of branching depends here on the short-scale parameters /3 and 7. 
Therefore, it cannot be predicted from purely energetic considerations. 

To interpret our results for accelerating cracks, it is useful to first derive an analytic 
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expression for the speed of steady-state cracks close to the Griffith threshold. Energy balance 
for steady-state propagation implies that the stored elastic energy ahead of the crack in 
excess of twice the surface energy must be dissipated in the process zone. This implies that 

dE „i f r^fd(l)^ ^ 



where the second equality follows from Eqs. (j21) with no energy flux through the boundaries. 
In a coordinate system {x' = x — vt, y' = y) translating with the tip, dip/dt = —vd(j)/dx'. 
Thus the steady-state velocity of rectilinear crack propagation is 

V ^ 27(G/G, - 1) 



where I = J dx'dy' {d(j)/dx')^ is a dimensionless integral factor of order unity that depends 
on the proflle of for a stationary crack G ^ Gc- This prediction is consistent with the 
simulation results that show a steeper increase of velocity with load for smaller /3. 

The fact that the T — v plots for different accelerating cracks in Fig.|21do not superimpose 
on the steady-state curve clearly shows that the speed of these cracks is not uniquely de- 
termined by the energy flow rate to the tip. In our model, accelerating cracks require more 
energy per unit length of advance than steadily propagating cracks. This extra energy is 
reversibly stored in the process zone to be either consumed by branching or radiated away 
later. This effect is exacerbated for large loads and causes the crack to decelerate before 
branching. This deceleration is marked by the overshoot of the velocity plotted as a function 
of tip position (Fig. H)) and the fact that the corresponding T — v plot becomes double- valued 
(Fig. 12)). A qualitatively similar overshoot of the velocity was observed experimentally for 
very large accelerations in brittle fracture of glass 3]; albeit for much larger systems than 
in the present simulations. 

To estimate when crack acceleration should be important in our model, we compute the 
change of the velocity over a distance comparable to the process zone. Crack's acceleration 
should be irrelevant when this change is small compared to the wave speed, or C, d{v/c)/dx <^ 
1. For cracks accelerating from rest, let us assume that Eq. gives the instantaneous speed 
of the accelerating crack if G is taken to mean the instantaneous elastic energy release rate 
F which varies as dV/dx ~ {G — Gc)/W as the crack extends into the stressed strip. This 
implies that the acceleration is small when 

W G-Gc 



We find that acceleration remains important in our simulations even for W/^ = 100 when 
f3 is of order unity. Since in experiments G can be up to a hundred times larger than 
Gc and (3 = ct/^ can be much smaller than unity, acceleration may be important even 
in macroscopic strips several orders of magnitude larger than the process zone size. We 
conclude that the validity of the continuum theory of brittle fracture depends not only on 
the ratio of the system size to the process zone size, but also on the importance of inertia 
relative to bond breaking dissipation measured here by (3. 

This conclusion does not conflict with the fact that this theory was recently validated 
by lattice simulations in the limit of zero dissipation In these simulations, the velocity 
jumps discontinuously to a fraction of the wave speed and the slope of the velocity-load 
relation is finite for velocities larger than this fraction. Therefore, 1/(3 should be replaced by 
d{v/c)/d{G/Gc) ~ 1 in the above estimate leading to Eq. ((7j). In contrast, for the velocity 
range from zero to a fraction of the wave speed pertinent to experiments, the continuum 
theory will break down for accelerating cracks in the /3 — limit of the present model for 
any finite system size. The crack dynamics is independent of (3 and well-defined in this limit. 

Clearly, the incorporation of more realistic microscopic mechanisms of failure in the 
present phase-field approach and the extension to mode I remain needed to make contact 
quantitatively with experiments. 
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